# Ryan Brutger and Alexandra Guisinger June 1, 2024
# Replicatin code for "Framing Layoffs" 
# Analysis run on macOS Catalina, MacBook Pro, 2.9 GHz Intel Core i5
# using R version 4.0.1 (2020-06-06) -- "See Things Now"

# The following code is for the survey experiment analysis

# ----------------------------------------------------------------------
# load all relevant themes and packages
# ----------------------------------------------------------------------

library("tidyverse")
library("ggplot2")
library("stargazer")
library("xtable")
library("texreg")
library("reshape2")
library("RItools")
library("estimatr")
library("ggpubr")
library("foreign")
library("dotwhisker")

# The following are the variable names for the experimental treatments:
# treat_m = Market Conditions
# treat_tb = Tariff Costs
# treat_tg = Tariff Competition
# treat_p = Pandemic

#Set your working directory
# setwd() 

#Load data
data <- read_csv("Replication_Data_Brutger_Guisinger_Political_Behavior-5-30-2024.csv")

# code gender identifier variables
data$women <- ifelse(data$gender==2, 1, 0)
data$men <- ifelse(data$gender==1, 1, 0)

#Baseline blame attribution in control, reported in note to Figure 3
table(data$blame_gov[data$AutoTreat==1])
121/(121+539) #0.183
table(data$blame_gm[data$AutoTreat==1])
334/(334+326) #0.506
table(data$blame_other[data$AutoTreat==1])
205/(205+455) #0.311

base.names <- c("Government", "General Motors", "Other")
base.opin <- c(.183, .506, .311)
color <- c("gray87","gray61", "gray29")

#Figure 3 of paper
barplot(base.opin, main="",  ylim=c(0,.6), 
        ylab="Proportion of Respondents", names.arg= base.names, col=color)
# Export 5x5

# For blame attribution reported in the text:
# blame attribution in control in US
table(data$blame_gov[data$US_study==1 & data$AutoTreat==1])
67/(67+249) #0.212
table(data$blame_gm[data$US_study==1 & data$AutoTreat==1])
144/(144+172) #0.456
table(data$blame_other[data$US_study==1 & data$AutoTreat==1])
105/(105+211) #0.332

# blame attribution in control in Canada
table(data$blame_gov[data$US_study==0 & data$AutoTreat==1])
54/(54+290) #0.157
table(data$blame_gm[data$US_study==0 & data$AutoTreat==1])
190/(190+154) #0.552
table(data$blame_other[data$US_study==0 & data$AutoTreat==1])
100/(100+244) #0.291

# Treatment effects on blame attribution used to calculate Figure 4 of paper
mod.gov <- lm(blame_gov ~ treat_m + treat_tb + treat_tg + treat_p , data=data)
mod.gm <- lm(blame_gm ~ treat_m + treat_tb + treat_tg + treat_p , data=data)
mod.other <- lm(blame_other ~ treat_m + treat_tb + treat_tg + treat_p , data=data)

#Figure 4 of paper
# Dot-and-whisker plots for binary variable
dw_blame <- dwplot(list(mod.gov, mod.gm, mod.other),
                   vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2,), # plot line at zero _behind_ coefs
                   dot_args = list(size=2)) %>%
  relabel_predictors(c(treat_m = "Market Conditions",                       
                       treat_tb = "Tariff-Costs", 
                       treat_tg = "Tariff-Competition",
                       treat_p = "Pandemic"
  ))+
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  theme_minimal() +
  theme(legend.position="bottom")+
  scale_colour_grey(name = "Dependent Variable",
                    labels = c("Blame Gov.", "Blame GM", "Blame Other"),                
    start = .3,
    end = .7)
print(dw_blame)
#ggsave(filename="dw_blame.pdf", plot=dw_blame, device="pdf", width=6, height=4, units="in")

# Following code is for Figure 5 of paper
# Treatment effects on blame attribution for those following the news

# Effects among those (not) following the auto news story
data.follow <- subset(data[data$follow.news==1 & is.na(data$follow.news)==FALSE, ])
data.not.follow <- subset(data[data$follow.news==0 & is.na(data$follow.news)==FALSE, ])
summary(data$follow.news) # 37% following

mod.gov.f <- lm(blame_gov ~ treat_m + treat_tb + treat_tg + treat_p , data=data.follow )
mod.gm.f <- lm(blame_gm ~ treat_m + treat_tb + treat_tg + treat_p , data=data.follow )
mod.other.f <- lm(blame_other ~ treat_m + treat_tb + treat_tg + treat_p , data=data.follow )

# Dot-and-whisker plot for panel (a) of Figure 5
dw_blame.f <- dwplot(list(mod.gov.f, mod.gm.f, mod.other.f),
                     vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2,), # plot line at zero _behind_ coefs
                     dot_args = list(size=2)) %>%
  relabel_predictors(c(treat_m = "Market Conditions",                       
                       treat_tb = "Tariff-Costs", 
                       treat_tg = "Tariff-Competition",
                       treat_p = "Pandemic"
  ))+
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  theme_minimal() +
  scale_color_hue(name = "Model",
                  labels = c("Blame Gov.", "Blame GM", "Blame Other"))+ 
  theme(legend.position="bottom")+
  scale_colour_grey(name = "Dependent Variable",
                    labels = c("Blame Gov.", "Blame GM", "Blame Other"),                
                    start = .3,
                    end = .7)
print(dw_blame.f)
#ggsave(filename="dw_blame_f.pdf", plot=dw_blame.f, device="pdf", width=6, height=3, units="in")

# Effects among those NOT following the auto news story
# Treatment effects on blame attribution NOT following the auto news story
mod.gov.nf <- lm(blame_gov ~ treat_m + treat_tb + treat_tg + treat_p , data=data.not.follow )
mod.gm.nf <- lm(blame_gm ~ treat_m + treat_tb + treat_tg + treat_p , data=data.not.follow )
mod.other.nf <- lm(blame_other ~ treat_m + treat_tb + treat_tg + treat_p , data=data.not.follow )

# Dot-and-whisker plot for panel (b) of Figure 5
dw_blame.nf <- dwplot(list(mod.gov.nf, mod.gm.nf, mod.other.nf),
                      vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2,), # plot line at zero _behind_ coefs
                      dot_args = list(size=2)) %>%
  relabel_predictors(c(treat_m = "Market Conditions",                       
                       treat_tb = "Tariff-Costs", 
                       treat_tg = "Tariff-Competition",
                       treat_p = "Pandemic"
  ))+
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  theme_minimal() +
  scale_color_hue(name = "Dependent Variable",
                  labels = c("Blame Gov.", "Blame GM", "Blame Other"))+ 
  theme(legend.position="bottom")+
  scale_colour_grey(name = "Dependent Variable",
                    labels = c("Blame Gov.", "Blame GM", "Blame Other"),                
                    start = .3,
                    end = .7)
print(dw_blame.nf)
#ggsave(filename="dw_blame_nf.pdf", plot=dw_blame.nf, device="pdf", width=6, height=3, units="in")

# Treatment effects on trade support
# create dichotomous trade support variable
data$auto_tr2<-ifelse(data$auto_tr>3, 1, 0)

mod.tr.supp <- lm(auto_tr ~ treat_m + treat_tb + treat_tg + treat_p , data=data)
summary(mod.tr.supp )# Tariff-Competition -0.09, p=0.05; Tariff-Costs 0.08, p = 0.08

# Trade support model with controls
mod.tr.supp.C <- lm(auto_tr ~ treat_m + treat_tb + treat_tg + treat_p + women + lib_conserv + education + age + US_study +income, data=data)

# For substantive percentage point changes, we use the dichotmous support variable
mod.tr.supp2 <- lm(auto_tr2 ~ treat_m + treat_tb + treat_tg + treat_p , data=data)
summary(mod.tr.supp2)# Tariff-Competition -3.1%, Tariff-Costs 5.1%, p=0.02

# Plot with sparse model and controls model - Figure 6 of the manuscript
dw_trade2 <- dwplot(list(mod.tr.supp, mod.tr.supp.C),
                    vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2,), # plot line at zero _behind_ coefs
                    dot_args = list(size=2)) %>%
  relabel_predictors(c(treat_m = "Market Conditions",                       
                       treat_tb = "Tariff-Costs", 
                       treat_tg = "Tariff-Competition",
                       treat_p = "Pandemic",
                       women = "Women",
                       lib_conserv = "Ideology",
                       education = "College Degree",
                       age = "Age",
                       US_study = "US Study",
                       income = "Income"
  ))+
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  theme_minimal() + 
  xlim(-.3,.3) +
  scale_color_hue(name = "Models:",
                  labels = c("Sparse Model", "Includes Controls"))+ 
  theme(legend.position="bottom")+
  scale_colour_grey(name = "Models",
                    labels = c("Sparse Model", "Includes Controls"),                
                    start = .3,
                    end = .7)
print(dw_trade2)
#ggsave(filename="dw_trade2.pdf", plot=dw_trade2, device="pdf", width=6, height=4, units="in")

##
##
## Code for Appendix
##
##

# Demographics reported in Table 1 of appendix

# DEMOGRAPHICS----
data$age_18to24 <- ifelse(data$age >=18 & data$age <25,1,0)
data$age_25to39 <- ifelse(data$age >=25 & data$age <40, 1,0)
data$age_40to59 <- ifelse(data$age >= 40 & data$age <60, 1,0)
data$age_60plus <- ifelse(data$age >=60, 1,0)

data$less50k <- ifelse(data$income<=2,1,0)
data$btw50_100k <- ifelse(data$income>2 & data$income<=4, 1,0)
data$btw100_150k <- ifelse(data$income>4 & data$income <=6, 1,0)
data$above150k <- ifelse(data$income>=7,1,0)

demo_names <- c("Age 18 to 24",
                "Age 25 to 39",
                "Age 40 to 59",
                "Age >50",
                "Female",
                "Household income $0 to $50,000",
                "Household income $50,001 to $100,000",
                "Household income $100,001 to $150,000",
                "Household income >$150,000",
                "Attended college")
demo_vars <- list(data$age_18to24,
                  data$age_25to39,
                  data$age_40to59,
                  data$age_60plus,
                  data$gender, # 1 = male, 2 = female
                  data$less50k,
                  data$btw50_100k,
                  data$btw100_150k,
                  data$above150k,
                  data$education)

data.US <- subset(data[data$US_study==1, ])
demo_vars_US <- list(data.US$age_18to24,
                     data.US$age_25to39,
                     data.US$age_40to59,
                     data.US$age_60plus,
                     data.US$gender, # 1 = male, 2 = female
                     data.US$less50k,
                     data.US$btw50_100k,
                     data.US$btw100_150k,
                     data.US$above150k,
                     data.US$education)

data.CA <- subset(data[data$US_study==0, ])
demo_vars_CA <- list(data.CA$age_18to24,
                     data.CA$age_25to39,
                     data.CA$age_40to59,
                     data.CA$age_60plus,
                     data.CA$gender, # 1 = male, 2 = female
                     data.CA$less50k,
                     data.CA$btw50_100k,
                     data.CA$btw100_150k,
                     data.CA$above150k,
                     data.CA$education)

pop_pct <- sprintf("%.3f",
                   round(c(# age breakdown from US Census Bureau 2019 Table 1: https://www.census.gov/data/tables/2019/demo/age-and-sex/2019-age-sex-composition.html (I ignore those <18)
                     0.1321,
                     0.2660,
                     0.3251,
                     0.2929,
                     .5097, # sex breakdown from US Census Bureau, 2019
                     .3707, #https://www.census.gov/data/tables/time-series/demo/income-poverty/cps-hinc/hinc-01.html
                     .2884,
                     .1555,
                     .1854,
                     0.6108 
                   ),digits=3)) # https://www.census.gov/content/census/en/data/tables/2019/demo/educational-attainment/cps-detailed-tables.html

prop.function <- function(data){
  prop <- as.numeric(prop.table(table(data))[2])
  prop <- sprintf("%.3f",round(prop,digits=3))
  return(prop)
}             

props_US <- lapply(demo_vars_US, prop.function)
props_CA <- lapply(demo_vars_CA, prop.function)

demo_props <- as.data.frame(cbind(demo_names, props_CA, props_US, pop_pct))
colnames(demo_props) <- c("Demographic", "Portion of Canadian Sample", "Portion of US Sample",  "U.S. Population")
print(demo_props, row.names=F)

print(xtable(demo_props,
             caption = "Study demographics. U.S. population information on age, sex, income, and education are from the Census Bureau and are for 2019. Partisan identification is from Pew and covers registered voters for 2018/9."),
      include.rownames=F)

# Sample sizes in each country: (manually added to table 1 of appendix)
length(data$US_study[data$US_study==1]) #2865 US sample
length(data$US_study[data$US_study==0]) #3118 Canadian Sample

# Table 2 of the appendix is manually generated from https://www150.statcan.gc.ca/

# Blame Attribution Regression Results
# Results reported in Table 3 of the appendix
stargazer(mod.gov, mod.gm, mod.other, omit.stat=c("LL", "ser", "f"), style="apsr", digits=3,
          title = "Blame Attribution", align = TRUE,
          covariate.labels = c("Market Conditions", "Tariff-Costs", "Tariff-Competition", "Pandemic"))

# Do Canadians and Americans respond differently?
# Interaction models reported in Table 4 of appendix
mod.gov.USint <- lm(blame_gov ~ treat_m*US_study + treat_tb*US_study + treat_tg*US_study + treat_p*US_study , data=data)
mod.gm.USint <- lm(blame_gm ~ treat_m*US_study  + treat_tb*US_study  + treat_tg*US_study  + treat_p*US_study  , data=data)
mod.other.USint <- lm(blame_other ~ treat_m*US_study  + treat_tb*US_study  + treat_tg*US_study  + treat_p*US_study , data=data)
stargazer(mod.gov.USint, mod.gm.USint, mod.other.USint) # order of rows in table are manually modified in appendix

# Effects among those (not) following the auto news story
data.follow <- subset(data[data$follow.news==1 & is.na(data$follow.news)==FALSE, ])
data.not.follow <- subset(data[data$follow.news==0 & is.na(data$follow.news)==FALSE, ])
summary(data$follow.news) # 37% following

# Do Canadians and Americans respond differently?
# Interaction models reported in Table 4 of appendix
mod.gov.USint <- lm(blame_gov ~ treat_m*US_study + treat_tb*US_study + treat_tg*US_study + treat_p*US_study , data=data)
mod.gm.USint <- lm(blame_gm ~ treat_m*US_study  + treat_tb*US_study  + treat_tg*US_study  + treat_p*US_study  , data=data)
mod.other.USint <- lm(blame_other ~ treat_m*US_study  + treat_tb*US_study  + treat_tg*US_study  + treat_p*US_study , data=data)
stargazer(mod.gov.USint, mod.gm.USint, mod.other.USint) # order of rows in table are manually modified in appendix

# Code prefered government policies
#Gov should provide job training & education support & tuition
data$train_edu <- data$should_train + data$should_educ
# Gov should provide unemployment benefits and wage wupplements
data$wage_ben <- data$should_unemp_ben + data$should_wage_supp

# Rescale so each measure takes a value from 0 to 1
rescale <- function(x){
  return((x-min(x,na.rm=TRUE))/(max(x-min(x,na.rm=TRUE),na.rm=TRUE)))
}
data$train_edu<- rescale(data$train_edu)
data$wage_ben<- rescale(data$wage_ben)

# Full sample interactions with US_Sample with controls reported in Table 5 of the appendix
mod.train.edu.C.Int <- lm(train_edu ~ treat_m*US_study + treat_tb*US_study + treat_tg*US_study + treat_p*US_study + women + lib_conserv + education + age + US_study +income, data=data)
mod.wage.ben.C.Int <- lm(wage_ben ~ treat_m*US_study + treat_tb*US_study + treat_tg*US_study + treat_p*US_study + women + lib_conserv + education + age + US_study +income , data=data)
stargazer(mod.train.edu.C.Int , mod.wage.ben.C.Int , omit.stat=c("LL", "ser", "f"), style="apsr", digits=3,
          title = "Heterogenous Effects on Support for Government Assistance", align = TRUE)

# Do those on the ideological right react differently?
# Interaction models reported in Table 6 of the appendix
mod.gov.IdeolInt <- lm(blame_gov ~ treat_m*ideology + treat_tb*ideology + treat_tg*ideology + treat_p*ideology , data=data)
mod.gm.IdeolInt <- lm(blame_gm ~ treat_m*ideology  + treat_tb*ideology  + treat_tg*ideology  + treat_p*ideology  , data=data)
mod.other.IdeolInt <- lm(blame_other ~ treat_m*ideology  + treat_tb*ideology  + treat_tg*ideology  + treat_p*ideology , data=data)
stargazer(mod.gov.IdeolInt, mod.gm.IdeolInt, mod.other.IdeolInt)

# Interaction model for US Respondents reported in Table 7 of appendix
mod.gov.IdeolIntUS <- lm(blame_gov ~ treat_m*ideology + treat_tb*ideology + treat_tg*ideology + treat_p*ideology , data=data[data$US_study==1, ])
mod.gm.IdeolIntUS <- lm(blame_gm ~ treat_m*ideology  + treat_tb*ideology  + treat_tg*ideology  + treat_p*ideology  , data=data[data$US_study==1, ])
mod.other.IdeolIntUS <- lm(blame_other ~ treat_m*ideology  + treat_tb*ideology  + treat_tg*ideology  + treat_p*ideology , data=data[data$US_study==1, ])
stargazer(mod.gov.IdeolIntUS, mod.gm.IdeolIntUS, mod.other.IdeolIntUS)

# Interaction model for Canadian Respondents reported in Table 8 of appendix
mod.gov.IdeolIntCan <- lm(blame_gov ~ treat_m*ideology + treat_tb*ideology + treat_tg*ideology + treat_p*ideology , data=data[data$US_study==0, ])
mod.gm.IdeolIntCan <- lm(blame_gm ~ treat_m*ideology  + treat_tb*ideology  + treat_tg*ideology  + treat_p*ideology  , data=data[data$US_study==0, ])
mod.other.IdeolIntCan <- lm(blame_other ~ treat_m*ideology  + treat_tb*ideology  + treat_tg*ideology  + treat_p*ideology , data=data[data$US_study==0, ])
stargazer(mod.gov.IdeolIntCan, mod.gm.IdeolIntCan, mod.other.IdeolIntCan)

#test the effect of ideology using a categorical variable
data$strongLib <- ifelse(data$ideology<3 ,1, 0)
data$moderate <- ifelse(data$ideology>2 & data$ideology<6,1, 0)
data$strongCon <- ifelse(data$ideology>5 ,1, 0)

# US Respondents Only - use categorical variables for stronger liberals, moderates, and stronger conservatives to see if they react differently
# Interaction models reported in Table 9 of the appendix
mod.gov.IdeolInt3US <- lm(blame_gov ~ treat_m*moderate + treat_tb*moderate + treat_tg*moderate + treat_p*moderate +treat_m*strongCon  + treat_tb*strongCon  + treat_tg*strongCon  + treat_p*strongCon  , data=data[data$US_study==1, ])
mod.gm.IdeolInt3US <- lm(blame_gm ~ treat_m*moderate  + treat_tb*moderate  + treat_tg*moderate  + treat_p*moderate  +treat_m*strongCon  + treat_tb*strongCon  + treat_tg*strongCon  + treat_p*strongCon  , data=data[data$US_study==1, ])
mod.other.IdeolInt3US <- lm(blame_other ~ treat_m*moderate  + treat_tb*moderate  + treat_tg*moderate  + treat_p*moderate +treat_m*strongCon  + treat_tb*strongCon  + treat_tg*strongCon  + treat_p*strongCon  , data=data[data$US_study==1, ])
stargazer(mod.gov.IdeolInt3US, mod.gm.IdeolInt3US, mod.other.IdeolInt3US)

# Canadaian Respondents Only - use categorical variables for stronger liberals, moderates, and stronger conservatives to see if they react differently
# Interaction models reported in Table 10 of the appendix
mod.gov.IdeolInt3Can <- lm(blame_gov ~ treat_m*moderate + treat_tb*moderate + treat_tg*moderate + treat_p*moderate +treat_m*strongCon  + treat_tb*strongCon  + treat_tg*strongCon  + treat_p*strongCon  , data=data[data$US_study==0, ])
mod.gm.IdeolInt3Can <- lm(blame_gm ~ treat_m*moderate  + treat_tb*moderate  + treat_tg*moderate  + treat_p*moderate  +treat_m*strongCon  + treat_tb*strongCon  + treat_tg*strongCon  + treat_p*strongCon  , data=data[data$US_study==0, ])
mod.other.IdeolInt3Can <- lm(blame_other ~ treat_m*moderate  + treat_tb*moderate  + treat_tg*moderate  + treat_p*moderate +treat_m*strongCon  + treat_tb*strongCon  + treat_tg*strongCon  + treat_p*strongCon  , data=data[data$US_study==0, ])
stargazer(mod.gov.IdeolInt3Can, mod.gm.IdeolInt3Can, mod.other.IdeolInt3Can)

# Testing for Hetergogenous Treatment Effects for Those (Not) Following the Auto Story
# Interaction models reported in Table 11 of the appendix
mod.gov.f.int <- lm(blame_gov ~ treat_m*follow.news + treat_tb*follow.news + treat_tg*follow.news + treat_p*follow.news, data=data)
mod.gm.f.int <- lm(blame_gm ~ treat_m*follow.news + treat_tb*follow.news + treat_tg*follow.news + treat_p*follow.news, data=data)
mod.other.f.int <- lm(blame_other ~ treat_m*follow.news + treat_tb*follow.news + treat_tg*follow.news + treat_p*follow.news, data=data)
stargazer(mod.gov.f.int, mod.gm.f.int, mod.other.f.int, omit.stat=c("LL", "ser", "f"), style="apsr", digits=3,
          title = "Blame Attribution", align = TRUE)

# Trade Support models reported in Table 12 of appendix
mod.tr.supp <- lm(auto_tr ~ treat_m + treat_tb + treat_tg + treat_p , data=data)
mod.tr.supp.C <- lm(auto_tr ~ treat_m + treat_tb + treat_tg + treat_p + women + lib_conserv + education + age + US_study +income, data=data)
stargazer(mod.tr.supp, mod.tr.supp.C,  covariate.labels = c("Market Conditions", "Tariff-Costs", "Tariff-Competition", "Pandemic", "Women", "Ideology", "College Degree", "Age", "US Study", "Income")) # For Table XXXX of the appendix

# Models on Support for Government Assistance Programs reported in Table 13 of appendix
mod.train.edu <- lm(train_edu ~ treat_m + treat_tb + treat_tg + treat_p , data=data)
mod.wage.ben <- lm(wage_ben ~ treat_m + treat_tb + treat_tg + treat_p , data=data)
stargazer(mod.train.edu,mod.wage.ben )

# Plot with sparse models for Figure 10 of appendix
dw_gov_pol <- dwplot(list(mod.train.edu, mod.wage.ben),
                     vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2,), # plot line at zero _behind_ coefs
                     dot_args = list(size=2)) %>%
  relabel_predictors(c(treat_m = "Market Conditions",                       
                       treat_tb = "Tariff-Costs", 
                       treat_tg = "Tariff-Competition",
                       treat_p = "Pandemic"
  ))+
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  theme_minimal() + 
  xlim(-.15, .15)+
  scale_color_hue(name = "Models:",
                  labels = c("Edu. & Training", "Wage Support"))+ 
  theme(legend.position="bottom")+
  scale_colour_grey(name = "Dependent Variable",
                    labels = c("Edu. & Training", "Wage Support"),                
                    start = .3,
                    end = .7)
print(dw_gov_pol)
#ggsave(filename="dw_gov_pol.pdf", plot=dw_gov_pol, device="pdf", width=6, height=4, units="in")

# Full sample with controls for Figure 11 of appendix
mod.train.edu.C <- lm(train_edu ~ treat_m + treat_tb + treat_tg + treat_p + women + lib_conserv + education + age + US_study +income, data=data)
mod.wage.ben.C <- lm(wage_ben ~ treat_m + treat_tb + treat_tg + treat_p + women + lib_conserv + education + age + US_study +income , data=data)

# Plot with controls for Figure 11 of appendix
dw_gov_pol2 <- dwplot(list(mod.train.edu.C, mod.wage.ben.C),
                      vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2,), # plot line at zero _behind_ coefs
                      dot_args = list(size=2)) %>%
  relabel_predictors(c(treat_m = "Market Conditions",                       
                       treat_tb = "Tariff-Costs", 
                       treat_tg = "Tariff-Competition",
                       treat_p = "Pandemic",
                       women = "Women",
                       lib_conserv = "Ideology",
                       education = "College Degree",
                       age = "Age",
                       US_study = "US Study",
                       income = "Income"
  ))+
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  theme_minimal() + 
  xlim(-.15, .15)+
  scale_color_hue(name = "Models:",
                  labels = c("Edu. & Training + Controls", "Wage Support + Controls"))+ 
  theme(legend.position="bottom")+
  scale_colour_grey(name = "Dependent Variable",
                    labels = c("Edu. & Training + Controls", "Wage Support + Controls"),                
                    start = .3,
                    end = .7)
print(dw_gov_pol2)
#ggsave(filename="dw_gov_pol2.pdf", plot=dw_gov_pol2, device="pdf", width=6, height=4, units="in")

# Testing for Heterogeneous Treaments E ects Based on Economic Vulnerability
# Interaction for those making less than 50k/year
data$less50k <- ifelse(data$income<=2,1,0)

# Models reported in Table 14 of appendix
mod.gov.lowInc <- lm(blame_gov ~ treat_m*less50k + treat_tb*less50k + treat_tg*less50k + treat_p*less50k , data=data)
mod.gm.lowInc <- lm(blame_gm ~ treat_m*less50k  + treat_tb*less50k  + treat_tg*less50k  + treat_p*less50k  , data=data)
mod.other.lowInc <- lm(blame_other ~ treat_m*less50k  + treat_tb*less50k  + treat_tg*less50k  + treat_p*less50k , data=data)
stargazer(mod.gov.lowInc, mod.gm.lowInc, mod.other.lowInc)

# If your employer was facing financial difficulties, would you be more or less likely to lose your job relative to other employees at the company?
# lower values mean higher likelihood of being laid off relative to peers
data$job_instability <- ifelse(data$Pre_Job_Stability<3, 1, 0)

# Models reported in Table 15 of the appendix
mod.gov.job.insec <- lm(blame_gov ~ treat_m*job_instability + treat_tb*job_instability + treat_tg*job_instability + treat_p*job_instability , data=data)
mod.gm.job.insec <- lm(blame_gm ~ treat_m*job_instability  + treat_tb*job_instability  + treat_tg*job_instability  + treat_p*job_instability  , data=data)
mod.other.job.insec <- lm(blame_other ~ treat_m*job_instability  + treat_tb*job_instability  + treat_tg*job_instability  + treat_p*job_instability , data=data)
stargazer(mod.gov.job.insec, mod.gm.job.insec, mod.other.job.insec)


# Do those with higher levels of education (college degree) respond differently? 
# Test education interactions with blame attribution
# Models reported in Table 16 of the appendix
mod.gov.EduInt <- lm(blame_gov ~ treat_m*education + treat_tb*education + treat_tg*education + treat_p*education , data=data)
mod.gm.EduInt <- lm(blame_gm ~ treat_m*education  + treat_tb*education  + treat_tg*education  + treat_p*education  , data=data)
mod.other.EduInt <- lm(blame_other ~ treat_m*education  + treat_tb*education  + treat_tg*education  + treat_p*education , data=data)
stargazer(mod.gov.EduInt, mod.gm.EduInt, mod.other.EduInt)

# Test education interactions with policy preferences
# Models reported in Table 17 of the appendix
mod.trade.EduInt <- lm(auto_tr ~ treat_m*education + treat_tb*education + treat_tg*education + treat_p*education , data=data)
mod.edu.EduInt <- lm(train_edu ~ treat_m*education  + treat_tb*education  + treat_tg*education  + treat_p*education  , data=data)
mod.wage.EduInt <- lm(wage_ben ~ treat_m*education  + treat_tb*education  + treat_tg*education  + treat_p*education , data=data)
stargazer(mod.trade.EduInt, mod.edu.EduInt, mod.wage.EduInt)

# Test ideology interactions with treatment for trade support and government assistance preferences
# models reported in Table 18 of the appendix
mod.tr.IdeolInt <- lm(auto_tr ~ treat_m*ideology + treat_tb*ideology + treat_tg*ideology + treat_p*ideology , data=data)
mod.edu.IdeolInt <- lm(train_edu ~ treat_m*ideology  + treat_tb*ideology  + treat_tg*ideology  + treat_p*ideology  , data=data)
mod.wage.IdeolInt <- lm(wage_ben ~ treat_m*ideology  + treat_tb*ideology  + treat_tg*ideology  + treat_p*ideology , data=data)
stargazer(mod.tr.IdeolInt, mod.edu.IdeolInt, mod.wage.IdeolInt)
